%pip install folium kneed prophet

import numpy as np                   
import pandas as pd   
import folium
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats          
import cvxpy as cp

from sklearn.metrics import mean_squared_error, r2_score
from folium.plugins import HeatMap
from pathlib import Path
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from pathlib import Path
from sklearn.cluster import KMeans, DBSCAN, AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import silhouette_score, calinski_harabasz_score, davies_bouldin_score
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.decomposition import PCA, KernelPCA
from sklearn.model_selection import GridSearchCV
from kneed import KneeLocator
from prophet import Prophet
from pandas import to_datetime
from prophet.diagnostics import cross_validation, performance_metrics

# Setting up inline plotting
%matplotlib inline

# Load Dataset
file_path = Path("") / "data" / "Cleaned_merged_SIPRI_Region_ACLED_starting2000.csv"     


# Load the CSV file
data = pd.read_csv(file_path)          

columns = data.info() # Create a Dataframe of Descriptive Statistics
columns             # Display the Dataframe of Descriptive Statistics

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 3716 entries, 0 to 3715
Data columns (total 29 columns):
 #   Column                              Non-Null Count  Dtype  
---  ------                              --------------  -----  
 0   Country                             3716 non-null   object 
 1   Year                                3716 non-null   int64  
 2   Expenditure-Share_of_Govt_spending  3533 non-null   float64
 3   Expenditure-Share_of_GDP            3602 non-null   float64
 4   Expenditure_Per_Capita              3604 non-null   float64
 5   Expenditure_Constant_2022           3673 non-null   float64
 6   Region                              3716 non-null   object 
 7   Subregion                           3716 non-null   object 
 8   index_level                         3716 non-null   int64  
 9   total_rank_from_avg_rank            3716 non-null   int64  
 10  avg_rank                            3716 non-null   float64
 11  Deadliness_rank                     3716 non-null   int64  
 12  Diffusion_rank                      3716 non-null   int64  
 13  Danger_rank                         3716 non-null   int64  
 14  Fragmentation_rank                  3716 non-null   int64  
 15  Deadliness_scaled_rank              3716 non-null   int64  
 16  Diffusion_scaled_rank               3716 non-null   int64  
 17  Danger_scaled_rank                  3716 non-null   int64  
 18  Fragmentation_scaled_rank           3716 non-null   int64  
 19  total_score_rank                    3716 non-null   int64  
 20  Deadliness_raw                      3716 non-null   int64  
 21  Diffusion_raw                       3716 non-null   float64
 22  Danger_raw                          3716 non-null   int64  
 23  Fragmentation_raw                   3716 non-null   int64  
 24  Deadliness_scaled                   3716 non-null   float64
 25  Diffusion_scaled                    3716 non-null   float64
 26  Danger_scaled                       3716 non-null   float64
 27  Fragmentation_scaled                3716 non-null   float64
 28  total_score                         3716 non-null   float64
dtypes: float64(11), int64(15), object(3)
memory usage: 842.0+ KB

# Descriptive Statistics
summary_stats = data.describe() # Create a Dataframe of Descriptive Statistics
summary_stats   # Display the Dataframe of Descriptive Statistics

# Drop the columns that are not relevant for the forecasting 
# model as mentioned by the user
columns_to_drop = [
    'Deadliness_scaled', 'Diffusion_scaled', 'Danger_scaled', 'Fragmentation_scaled',
    'Deadliness_rank', 'Diffusion_rank', 'Danger_rank', 'Fragmentation_rank'
]
data_cleaned = data.drop(columns=columns_to_drop, errors='ignore')

# Display the cleaned dataset to check the changes
data_cleaned.head()

# Check for missing values in the cleaned dataset
missing_values = data_cleaned.isna().sum()

# Display the columns with missing values and the count of missing values
missing_values[missing_values > 0]

Expenditure-Share_of_Govt_spending    183
Expenditure-Share_of_GDP              114
Expenditure_Per_Capita                112
Expenditure_Constant_2022              43
dtype: int64

# Create a copy of the dataframe to avoid modifying the original data
data_cleaned = data_cleaned.copy()

# Handling missing values in "Expenditure_Constant_2022" using forward and backward filling
data_cleaned['Expenditure_Constant_2022'] = data_cleaned['Expenditure_Constant_2022'].ffill().bfill()

# Handling missing values in other numerical columns using median imputation
for column in ['Expenditure-Share_of_Govt_spending', 'Expenditure-Share_of_GDP', 
               'Expenditure_Per_Capita', 'Expenditure_Constant_2022']:
    median_value = data_cleaned[column].median()
    data_cleaned[column] = data_cleaned[column].fillna(median_value)

# Verify that there are no missing values left in the dataset
missing_values_after_imputation = data_cleaned.isna().sum()
missing_values_after_imputation[missing_values_after_imputation > 0]

Series([], dtype: int64)

data_for_clustering = data.copy()

# Checking the percentage of missing values in each column
missing_values = data_for_clustering.isnull().mean() * 100
missing_values = missing_values[missing_values > 0].sort_values(ascending=False)
missing_values

Expenditure-Share_of_Govt_spending    4.924650
Expenditure-Share_of_GDP              3.067815
Expenditure_Per_Capita                3.013994
Expenditure_Constant_2022             1.157158
dtype: float64

# Impute missing values in the specified columns using the median
columns_to_impute = [
    "Expenditure-Share_of_Govt_spending", 
    "Expenditure-Share_of_GDP", 
    "Expenditure_Per_Capita", 
    "Expenditure_Constant_2022"
]

# Impute missing values in the specified columns using the median
for column in columns_to_impute:
    data_for_clustering[column] = data_for_clustering[column]\
        .fillna(data_for_clustering[column].median())

# Verify if there are any missing values left in these columns
remaining_missing = data_for_clustering[columns_to_impute].isnull().sum()
remaining_missing

Expenditure-Share_of_Govt_spending    0
Expenditure-Share_of_GDP              0
Expenditure_Per_Capita                0
Expenditure_Constant_2022             0
dtype: int64

# Creating a mapping for renaming the columns
column_rename_map = {
    "Expenditure-Share_of_Govt_spending": "GovtSpending",
    "Expenditure-Share_of_GDP": "GDP",
    "Expenditure_Per_Capita": "PerCapita",
    "Expenditure_Constant_2022": "ExpenditureConstant",
    "total_score": "TotalScore",
    "Deadliness_raw": "Deadliness",
    "Diffusion_raw": "Diffusion",
    "Danger_raw": "Danger",
    "Fragmentation_raw": "Fragmentation"
}

# Renaming the columns in the dataset
data_for_clustering.rename(columns=column_rename_map, inplace=True)

# Plotting histograms for numerical columns
numerical_columns = data.select_dtypes(include=['float64', 'int64']).columns
data[numerical_columns].hist(figsize=(15, 15), bins=20)
plt.suptitle('Histograms of Numerical Data', size=16)
plt.show()

# Plotting bar graphs for categorical columns
categorical_columns = data.select_dtypes(include=['object']).columns

plt.figure(figsize=(15, 20))
for i, column in enumerate(categorical_columns):
    plt.subplot(len(categorical_columns), 1, i + 1)
    sns.countplot(data=data, x=column)
    plt.xticks(rotation=45)
    plt.title(f'Bar Graph of {column}')

plt.tight_layout()
plt.show()

# Generating a correlation heatmap for numerical data
# Calculate the correlation matrix for numerical columns
numerical_columns_cleaned = data_cleaned.select_dtypes(include=['float64', 'int64']).columns
correlation_matrix = data_cleaned[numerical_columns_cleaned].corr()

# Plotting the heatmap using seaborn
plt.figure(figsize=(15, 10))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Heatmap of Numerical Data', size=16)
plt.show()

# Identifying highly correlated features (absolute correlation > 0.7)
correlation_matrix = data[numerical_columns].corr()
highly_correlated_pairs = correlation_matrix.abs().unstack().sort_values(ascending=False)
highly_correlated_pairs = highly_correlated_pairs[highly_correlated_pairs < 1.0]  
# Exclude self-correlation

# Get the top correlated pairs
top_correlated_pairs = highly_correlated_pairs[highly_correlated_pairs > 0.7].index

# Plotting scatter plots for highly correlated feature pairs
plt.figure(figsize=(15, 10))
plot_index = 1
for (feature1, feature2) in top_correlated_pairs:
    if plot_index > 6:  # Limit the number of scatter plots to display
        break
    plt.subplot(2, 3, plot_index)
    sns.scatterplot(data=data, x=feature1, y=feature2)
    plt.title(f'Scatter Plot: {feature1} vs {feature2}')
    plot_index += 1

plt.tight_layout()
plt.show()

# Updating the numerical_columns list with the new names
numerical_columns = list(column_rename_map.values())

# Calculate the correlation matrix
correlation_matrix = data_for_clustering[numerical_columns].corr()

# Plot the correlation heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f", cmap="coolwarm", \
            vmin=-1, vmax=1, linewidths=0.5)
plt.title("Correlation Matrix of Numerical Attributes", fontsize=16)
plt.tight_layout()
plt.show()

# Calculating the mean of numerical attributes grouped by 'Region' and 'Subregion'
region_patterns = data_for_clustering.groupby("Region")[numerical_columns].mean()
subregion_patterns = data_for_clustering.groupby("Subregion")[numerical_columns].mean()

# Plotting the region-level patterns as a heatmap
plt.figure(figsize=(12, 8))
sns.heatmap(region_patterns, annot=True, fmt=".2f", cmap="coolwarm", linewidths=0.5)
plt.title("Region-Level Patterns", fontsize=16)
plt.ylabel("Region", fontsize=12)
plt.xlabel("Attributes", fontsize=12)
plt.tight_layout()
plt.show()

# Recalculating conflict metric averages at the region and subregion levels for a clearer reassessment
conflict_metrics = ["TotalScore", "Deadliness", "Diffusion", "Danger", "Fragmentation"]

# Region-level conflict metrics mean
region_conflict_patterns = data_for_clustering.groupby("Region")[conflict_metrics].mean()

# Subregion-level conflict metrics mean
subregion_conflict_patterns = data_for_clustering.groupby("Subregion")[conflict_metrics].mean()

# Plotting the region-level conflict metrics as a heatmap
plt.figure(figsize=(12, 8))
sns.heatmap(region_conflict_patterns, annot=True, fmt=".2f", cmap="coolwarm", linewidths=0.5)
plt.title("Region-Level Conflict Metrics", fontsize=16)
plt.ylabel("Region", fontsize=12)
plt.xlabel("Conflict Metrics", fontsize=12)
plt.tight_layout()
plt.show()

# Ensure aggregation is applied only to numeric columns for countries
numeric_columns = data_for_clustering.select_dtypes(include=['number']).columns
country_summary = data_for_clustering.groupby('Country', as_index=False)[numeric_columns].mean()

# Define exact latitude and longitude for the given top countries
country_coordinates = {
    "United States of America": [39.8283, -98.5795],
    "China": [35.8617, 104.1954],
    "Russia": [61.5240, 105.3188],
    "United Kingdom": [55.3781, -3.4360],
    "Saudi Arabia": [23.8859, 45.0792],
    "India": [20.5937, 78.9629],
    "France": [46.6034, 1.8883],
    "Germany": [51.1657, 10.4515],
    "Japan": [36.2048, 138.2529],
    "Korea, South": [35.9078, 127.7669]
}

# Filter the dataset for the top countries
top_countries = data_for_clustering[data_for_clustering['Country'].isin(country_coordinates.keys())]

# Aggregate the numeric columns for the selected countries
numeric_columns = top_countries.select_dtypes(include=['number']).columns
top_countries_summary = top_countries.groupby('Country', as_index=False)[numeric_columns].mean()

# Add latitude and longitude for the selected countries
top_countries_summary['Latitude'] = top_countries_summary['Country']\
    .map(lambda x: country_coordinates[x][0])
top_countries_summary['Longitude'] = top_countries_summary['Country']\
    .map(lambda x: country_coordinates[x][1])

# Prepare data for the HeatMap layer
heat_data = top_countries_summary[['Latitude', 'Longitude', 'ExpenditureConstant']]\
    .dropna().values.tolist()

# Create a base map
m = folium.Map(location=[15, 0], zoom_start=2, png_enabled=True)

# Add HeatMap for the top countries
HeatMap(heat_data, radius=30, max_zoom=10, gradient={0.2: 'blue', 0.4: 'lime', \
                                                     0.6: 'yellow', 1: 'red'}).add_to(m)

# Add markers for each country
# for _, row in top_countries_summary.iterrows():
#     folium.Marker(
#         location=[row['Latitude'], row['Longitude']],
#         popup=f"Country: {row['Country']}<br>Expenditure: {row['ExpenditureConstant']:.2f}",
#         icon=folium.Icon(color="blue", icon="info-sign")
#     ).add_to(m)
    
m

# Selecting relevant features for forecasting
# We will primarily use 'Year' and 'Expenditure_Constant_2022' for the time-series analysis.
forecast_data = data_cleaned[['Country', 'Year', 'Expenditure_Constant_2022']]

# Sorting the data by Country and Year to prepare for time-series modeling
forecast_data = forecast_data.sort_values(by=['Country', 'Year'])

# Display the prepared data for forecasting
forecast_data.head()

# Split the data into training and testing sets for each country
# We will create a dictionary to hold the train-test split for each country
countries = forecast_data['Country'].unique()
train_test_data = {}

for country in countries:
    country_data = forecast_data[forecast_data['Country'] == country]
    X = country_data[['Year']]
    y = country_data['Expenditure_Constant_2022']
    
    # Split the data (80% training, 20% testing)
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)
    
    train_test_data[country] = {
        'X_train': X_train,
        'X_test': X_test,
        'y_train': y_train,
        'y_test': y_test
    }

# Display the train-test split information for a sample country
sample_country = countries[37]
print(sample_country)
train_test_data[sample_country]

Cuba

{'X_train':       Year
 1587  2003
 1588  2004
 1589  2005
 1590  2006
 1591  2007
 1592  2008
 1593  2009
 1594  2010
 1595  2011
 1596  2012
 1597  2013
 1598  2014,
 'X_test':       Year
 1599  2015
 1600  2016
 1601  2017
 1602  2018,
 'y_train': 1587    1290.519382
 1588    1290.519382
 1589    1290.519382
 1590    1290.519382
 1591    1290.519382
 1592    1290.519382
 1593    1290.519382
 1594    1290.519382
 1595    1290.519382
 1596    1290.519382
 1597    1290.519382
 1598    1290.519382
 Name: Expenditure_Constant_2022, dtype: float64,
 'y_test': 1599    1290.519382
 1600    1290.519382
 1601    1290.519382
 1602    1290.519382
 Name: Expenditure_Constant_2022, dtype: float64}

# Initialize a dictionary to store the models for each country
models = {}

# Train a Linear Regression model for each country
for country in countries:
    # Extract train-test data for the current country
    X_train = train_test_data[country]['X_train']
    y_train = train_test_data[country]['y_train']
    
    # Initialize and train the Linear Regression model
    model = LinearRegression()
    model.fit(X_train, y_train)
    
    # Store the trained model
    models[country] = model

# Forecast the next 5 and 10 years for each country using the trained models
forecast_years_5 = pd.DataFrame({'Year': np.arange(2022, 2027)})
forecast_years_10 = pd.DataFrame({'Year': np.arange(2022, 2032)})

# Create a dictionary to store forecasts for each country
forecasts = {}
for country in countries:
    # Use the trained model to make predictions for the next 5 and 10 years
    model = models[country]
    forecast_5 = model.predict(forecast_years_5)
    forecast_10 = model.predict(forecast_years_10)
  
    forecasts[country] = {
        '5_year_forecast': forecast_5,
        '10_year_forecast': forecast_10
    }
# Display the forecast for the first country as a sample
forecasts[sample_country]

{'5_year_forecast': array([1290.519382, 1290.519382, 1290.519382, 1290.519382, 1290.519382]),
 '10_year_forecast': array([1290.519382, 1290.519382, 1290.519382, 1290.519382, 1290.519382,
        1290.519382, 1290.519382, 1290.519382, 1290.519382, 1290.519382])}

# Get the most recent "Expenditure_Constant_2022" for each country
latest_expenditure = forecast_data.sort_values(['Country', 'Year'])\
    .groupby('Country').last()['Expenditure_Constant_2022']

# Sort and select the top 15 countries
top_15_countries = latest_expenditure.sort_values(ascending=False).head(15).index

# Plotting the 5-year forecast for the top 15 countries
plt.figure(figsize=(15, 10))
for country in top_15_countries:
    years = np.arange(2022, 2027)
    expenditure_forecast = forecasts[country]['5_year_forecast']
    plt.plot(years, expenditure_forecast, label=country)

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('5-Year Forecast for the Top 15 Countries by Expenditure_Constant_2022')
plt.legend(loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Plotting the 10-year forecast for the top 15 countries
plt.figure(figsize=(15, 10))
for country in top_15_countries:
    years = np.arange(2022, 2032)
    expenditure_forecast = forecasts[country]['10_year_forecast']
    plt.plot(years, expenditure_forecast, label=country)

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('10-Year Forecast for the Top 15 Countries by Expenditure_Constant_2022')
plt.legend(loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Sort the top 15 countries by the most recent expenditure level for legend ranking
latest_expenditure_sorted = latest_expenditure[top_15_countries].sort_values(ascending=False)

plt.figure(figsize=(15, 10))

country_lines = []

for country in latest_expenditure_sorted.index:
    # Get the historical data for the last 10 years for the selected country
    country_data = forecast_data[forecast_data['Country'] == country]
    historical_years = country_data['Year'].tail(10)
    historical_expenditure = country_data['Expenditure_Constant_2022'].tail(10)

    # Plot the historical data with a solid line
    historical_line, = plt.plot(
        historical_years, historical_expenditure, linestyle='-', label=f'{country} (History)'
    )
    
    # Plot the forecast for the next 10 years with a dashed line, using the same color
    plt.plot(
        np.arange(2022, 2032), forecasts[country]['10_year_forecast'], 
        linestyle='--', color=historical_line.get_color(), label=f'{country} (Forecast)'
    )
    
    # Store reference for legend sorting
    country_lines.append((country, historical_line))

# Sort legend entries by the latest expenditure level
sorted_legend_handles = [line for country, line in \
                         sorted(country_lines, key=lambda x: latest_expenditure_sorted[country], reverse=True)]

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('Historical and 10-Year Forecast for the Top 15 Countries by Expenditure_Constant_2022')

# Placing the legend outside to the left of the plot and ranked by most recent expenditure level
plt.legend(handles=sorted_legend_handles, loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Displaying the formulas (coefficients and intercepts) of the Linear Regression models for each country

# Store the forecast formulas for each of the top 15 countries
forecast_formulas = {}

for country in top_15_countries:
    # Extract the model for the current country
    model = models[country]
    
    # Get the coefficient (slope) and intercept of the linear regression model
    coefficient = model.coef_[0]
    intercept = model.intercept_
    
    # Store the formula
    forecast_formulas[country] = f"Expenditure_Constant_2022 = {coefficient:.2f} * Year + {intercept:.2f}"

# Display the forecast formulas for the top 15 countries
forecast_formulas

{'United States of America': 'Expenditure_Constant_2022 = 11981.00 * Year + -23272546.56',
 'China': 'Expenditure_Constant_2022 = 11618.69 * Year + -23207588.84',
 'Russia': 'Expenditure_Constant_2022 = 3464.69 * Year + -6903680.88',
 'Somalia': 'Expenditure_Constant_2022 = 0.00 * Year + 126473.35',
 'India': 'Expenditure_Constant_2022 = 2265.94 * Year + -4503217.61',
 'Saudi Arabia': 'Expenditure_Constant_2022 = 3301.45 * Year + -6574677.22',
 'United Kingdom': 'Expenditure_Constant_2022 = 81.02 * Year + -100966.65',
 'Ukraine': 'Expenditure_Constant_2022 = 194.25 * Year + -386780.27',
 'Germany': 'Expenditure_Constant_2022 = 135.10 * Year + -229344.53',
 'France': 'Expenditure_Constant_2022 = 174.92 * Year + -303368.71',
 'Japan': 'Expenditure_Constant_2022 = 43.12 * Year + -46864.64',
 'Korea, South': 'Expenditure_Constant_2022 = 1086.85 * Year + -2153701.90',
 'Italy': 'Expenditure_Constant_2022 = -558.10 * Year + 1151893.60',
 'Australia': 'Expenditure_Constant_2022 = 719.24 * Year + -1423274.48',
 'Israel': 'Expenditure_Constant_2022 = 454.85 * Year + -896067.90'}

# Initialize a dictionary to store the Random Forest models for each country
rf_models = {}

# Train a Random Forest Regressor for each country
for country in top_15_countries:
    # Extract train-test data for the current country
    X_train = train_test_data[country]['X_train']
    y_train = train_test_data[country]['y_train']
    
    # Initialize and train the Random Forest Regressor
    rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
    rf_model.fit(X_train, y_train)
    
    # Store the trained model
    rf_models[country] = rf_model

# Forecast the next 10 years for each country using the trained Random Forest models
rf_forecasts = {}

for country in top_15_countries:
    # Use the trained model to make predictions for the next 10 years
    rf_model = rf_models[country]
    forecast_10 = rf_model.predict(forecast_years_10)
    
    rf_forecasts[country] = {
        '10_year_forecast': forecast_10
    }

# Display the forecast for one of the top 15 countries
sample_country = top_15_countries[0]
rf_forecasts[sample_country]

{'10_year_forecast': array([789481.619213, 789481.619213, 789481.619213, 789481.619213,
        789481.619213, 789481.619213, 789481.619213, 789481.619213,
        789481.619213, 789481.619213])}

# Plotting the historical data (last 10 years) and the 10-year Random Forest 
# forecast for the top 15 countries with consistent colors

plt.figure(figsize=(15, 10))

for country in sorted(top_15_countries):
    # Get the historical data for the last 10 years for the selected country
    country_data = forecast_data[forecast_data['Country'] == country]
    historical_years = country_data['Year'].tail(10)
    historical_expenditure = country_data['Expenditure_Constant_2022'].tail(10)

    # Plot the historical data with a solid line
    historical_line, = plt.plot(
        historical_years, historical_expenditure, linestyle='-', label=f'{country} (History)'
    )
    
    # Plot the Random Forest forecast for the next 10 years with a dashed line, using the same color
    forecast_years = np.arange(2022, 2032)
    expenditure_forecast = rf_forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, expenditure_forecast, linestyle='--', \
            color=historical_line.get_color(), label=f'{country} (RF Forecast)'
    )

# Sort legend entries by the latest expenditure level
sorted_legend_handles = [line for country, line in \
                         sorted(country_lines, key=lambda x: latest_expenditure_sorted[x[0]], reverse=True)]

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('Historical and 10-Year Random Forest Forecast for the Top 15 Countries by Expenditure_Constant_2022')

# Placing the legend outside to the left of the plot and ranked by most recent expenditure level
plt.legend(handles=sorted_legend_handles, loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Plotting the historical data, Linear Regression forecast, and Random Forest 
# forecast for the top 15 countries with consistent colors

plt.figure(figsize=(15, 10))

for country in sorted(top_15_countries):
    # Get the historical data for the last 10 years for the selected country
    country_data = forecast_data[forecast_data['Country'] == country]
    historical_years = country_data['Year'].tail(10)
    historical_expenditure = country_data['Expenditure_Constant_2022'].tail(10)

    # Plot the historical data with a solid line
    historical_line, = plt.plot(
        historical_years, historical_expenditure, linestyle='-', label=f'{country} (History)'
    )
    
    # Plot the Linear Regression forecast for the next 10 years with a dashed line
    linear_forecast_years = np.arange(2022, 2032)
    linear_expenditure_forecast = forecasts[country]['10_year_forecast']
    plt.plot(
        linear_forecast_years, linear_expenditure_forecast, linestyle='--', color=historical_line.get_color(), 
        label=f'{country} (Linear Forecast)'
    )
    
    # Plot the Random Forest forecast for the next 10 years with a dotted line
    rf_expenditure_forecast = rf_forecasts[country]['10_year_forecast']
    plt.plot(
        linear_forecast_years, rf_expenditure_forecast, linestyle=':', color=historical_line.get_color(), 
        label=f'{country} (RF Forecast)'
    )

# Sort legend entries by the latest expenditure level
sorted_legend_handles = [line for country, line in \
                         sorted(country_lines, key=lambda x: latest_expenditure_sorted[x[0]], reverse=True)]

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('Historical, Linear Regression Forecast, and Random Forest Forecast for the Top 15 Countries')

# Placing the legend outside to the left of the plot and ranked by most recent expenditure level
plt.legend(handles=sorted_legend_handles, loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Prepare a list to store the summary of the Random Forest model's intercept 
# (using the mean), coefficient approximation (feature importance), and R^2 for each country
rf_model_summary = []

# Calculate the intercept (mean value), approximate coefficient, and R^2 score for each of the top 15 countries
for country in top_15_countries:
    # Extract the Random Forest model and train-test data for the current country
    rf_model = rf_models[country]
    X_train = train_test_data[country]['X_train']
    y_train = train_test_data[country]['y_train']
    X_test = train_test_data[country]['X_test']
    y_test = train_test_data[country]['y_test']
    
    # Calculate the intercept as the mean of the target value
    intercept = y_train.mean()
    
    # Approximate the "coefficient" as feature importance (since Random Forest doesn't directly provide coefficients)
    feature_importance = rf_model.feature_importances_[0]
    
    # Calculate the R^2 score on the test data
    y_pred = rf_model.predict(X_test)
    r_squared = r2_score(y_test, y_pred)
    
    # Append the summary to the list
    rf_model_summary.append({
        'Country': country,
        'Approximate Coefficient (Feature Importance)': feature_importance,
        'Intercept (Mean Value)': intercept,
        'R^2 Score': r_squared
    })

# Create a DataFrame to display the model summary
rf_model_summary_df = pd.DataFrame(rf_model_summary)

rf_model_summary_df

# Prepare a list to store the comparison of R^2 scores for Linear Regression and Random Forest for each country
model_comparison = []

# Calculate the R^2 score for each of the top 15 countries for both Linear Regression and Random Forest models
for country in top_15_countries:
    # Extract train-test data for the current country
    X_test = train_test_data[country]['X_test']
    y_test = train_test_data[country]['y_test']
    
    # Linear Regression R^2 Score
    linear_model = models[country]
    y_pred_linear = linear_model.predict(X_test)
    r_squared_linear = r2_score(y_test, y_pred_linear)
    
    # Random Forest R^2 Score
    rf_model = rf_models[country]
    y_pred_rf = rf_model.predict(X_test)
    r_squared_rf = r2_score(y_test, y_pred_rf)
    
    # Append the comparison to the list
    model_comparison.append({
        'Country': country,
        'Linear R^2 Score': r_squared_linear,
        'Random Forest R^2 Score': r_squared_rf
    })

# Create a DataFrame to display the model comparison
model_comparison_df = pd.DataFrame(model_comparison)

model_comparison_df

# Initialize a dictionary to store the Gradient Boosting models for each country
gbr_models = {}

# Train a Gradient Boosting Regressor for each country
for country in top_15_countries:
    # Extract train-test data for the current country
    X_train = train_test_data[country]['X_train']
    y_train = train_test_data[country]['y_train']
    
    # Initialize and train the Gradient Boosting Regressor
    gbr_model = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, max_depth=3, random_state=42)
    gbr_model.fit(X_train, y_train)
    
    # Store the trained model
    gbr_models[country] = gbr_model

# Forecast the next 10 years for each country using the trained Gradient Boosting models
gbr_forecasts = {}

for country in top_15_countries:
    # Use the trained model to make predictions for the next 10 years
    gbr_model = gbr_models[country]
    forecast_10 = gbr_model.predict(forecast_years_10)
    
    gbr_forecasts[country] = {
        '10_year_forecast': forecast_10
    }

# Display the forecast for the first country as a sample
gbr_forecasts[sample_country]

{'10_year_forecast': array([795404.95804064, 795404.95804064, 795404.95804064, 795404.95804064,
        795404.95804064, 795404.95804064, 795404.95804064, 795404.95804064,
        795404.95804064, 795404.95804064])}

from sklearn.ensemble import GradientBoostingRegressor

# Initialize a dictionary to store the Gradient Boosting models for each country
gbr_models = {}

# Train a Gradient Boosting Regressor for each country
for country in top_15_countries:
    # Extract train-test data for the current country
    X_train = train_test_data[country]['X_train']
    y_train = train_test_data[country]['y_train']
    
    # Initialize and train the Gradient Boosting Regressor
    gbr_model = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, max_depth=3, random_state=42)
    gbr_model.fit(X_train, y_train)
    
    # Store the trained model
    gbr_models[country] = gbr_model

# Forecast the next 10 years for each country using the trained Gradient Boosting models
gbr_forecasts = {}

for country in top_15_countries:
    # Use the trained model to make predictions for the next 10 years
    gbr_model = gbr_models[country]
    forecast_10 = gbr_model.predict(forecast_years_10)
    
    gbr_forecasts[country] = {
        '10_year_forecast': forecast_10
    }

# Display the forecast for the first country as a sample
gbr_forecasts[sample_country]

{'10_year_forecast': array([795404.95804064, 795404.95804064, 795404.95804064, 795404.95804064,
        795404.95804064, 795404.95804064, 795404.95804064, 795404.95804064,
        795404.95804064, 795404.95804064])}

# Plotting the historical data, Linear Regression forecast, Random Forest forecast, 
# and Gradient Boosting Regressor forecast for the top 15 countries

plt.figure(figsize=(15, 10))

for country in sorted(top_15_countries):
    # Get the historical data for the last 10 years for the selected country
    country_data = forecast_data[forecast_data['Country'] == country]
    historical_years = country_data['Year'].tail(10)
    historical_expenditure = country_data['Expenditure_Constant_2022'].tail(10)

    # Plot the historical data with a solid line
    historical_line, = plt.plot(
        historical_years, historical_expenditure, linestyle='-', label=f'{country} (History)'
    )
    
    # Plot the Linear Regression forecast for the next 10 years with a dashed line
    forecast_years = np.arange(2022, 2032)
    linear_expenditure_forecast = forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, linear_expenditure_forecast, linestyle='--', color=historical_line.get_color(), 
        label=f'{country} (Linear Forecast)'
    )
    
    # Plot the Random Forest forecast for the next 10 years with a dotted line
    rf_expenditure_forecast = rf_forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, rf_expenditure_forecast, linestyle=':', color=historical_line.get_color(), 
        label=f'{country} (RF Forecast)'
    )
    
    # Plot the Gradient Boosting Regressor forecast for the next 10 years with a dash-dot line
    gbr_expenditure_forecast = gbr_forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, gbr_expenditure_forecast, linestyle='-.', color=historical_line.get_color(), 
        label=f'{country} (GBR Forecast)'
    )

# Sort legend entries by the latest expenditure level
sorted_legend_handles = [line for country, line in \
                         sorted(country_lines, key=lambda x: latest_expenditure_sorted[x[0]], reverse=True)]

plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.title('Historical, Linear Regression, Random Forest, and Gradient Boosting Forecast for the Top 15 Countries')

# Placing the legend outside to the left of the plot and ranked by most recent expenditure level
plt.legend(handles=sorted_legend_handles, loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which="both", linestyle='--')
plt.show()

# Plotting the historical data, Linear Regression forecast, Random Forest 
# forecast, and Gradient Boosting Regressor forecast for the top 15 countries

plt.figure(figsize=(15, 10))

country_lines = []  # To store line references for sorting later

for country in sorted(top_15_countries):
    # Get the historical data for the last 10 years for the selected country
    country_data = forecast_data[forecast_data['Country'] == country]
    historical_years = country_data['Year'].tail(10)
    historical_expenditure = country_data['Expenditure_Constant_2022'].tail(10)

    # Plot the historical data with a solid line
    historical_line, = plt.plot(
        historical_years, historical_expenditure, linestyle='-', label=f'{country} (History)'
    )
    country_lines.append((country, historical_line))  # Storing line reference for legend sorting

    # Plot the Linear Regression forecast for the next 10 years with a dashed line
    forecast_years = np.arange(2022, 2032)
    linear_expenditure_forecast = forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, linear_expenditure_forecast, linestyle='--', color=historical_line.get_color(), 
        label=f'{country} (Linear Forecast)'
    )
    
    # Plot the Random Forest forecast for the next 10 years with a dotted line
    rf_expenditure_forecast = rf_forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, rf_expenditure_forecast, linestyle=':', color=historical_line.get_color(), 
        label=f'{country} (RF Forecast)'
    )
    
    # Plot the Gradient Boosting Regressor forecast for the next 10 years with a dash-dot line
    gbr_expenditure_forecast = gbr_forecasts[country]['10_year_forecast']
    plt.plot(
        forecast_years, gbr_expenditure_forecast, linestyle='-.', color=historical_line.get_color(), 
        label=f'{country} (GBR Forecast)'
    )

# Sort legend entries by the latest expenditure level
sorted_legend_handles = [line for country, line in \
                         sorted(country_lines, key=lambda x: historical_expenditure.iloc[-1], reverse=True)]

plt.xlabel('Year')
plt.ylabel('Expenditure (Constant 2022 Dollars)')
plt.yscale('log')  # Setting the y-axis to logarithmic scale
plt.title('Historical, Linear Regression, Random Forest, and Gradient Boosting Forecast for the Top 15 Countries')

# Placing the legend outside to the left of the plot and ranked by most recent expenditure level
plt.legend(handles=sorted_legend_handles, loc='upper left', bbox_to_anchor=(-0.3, 1), fontsize='small', ncol=1)
plt.grid(True, which='both', linestyle='--')  # Adding gridlines for better readability
plt.show()

# Prepare a list to store the summary of the Gradient Boosting model's R^2 for each country
gbr_model_summary = []

# Calculate the R^2 score for each of the top 15 countries for the Gradient Boosting model
for country in top_15_countries:
    # Extract train-test data for the current country
    X_test = train_test_data[country]['X_test']
    y_test = train_test_data[country]['y_test']
    
    # Gradient Boosting Regressor R^2 Score
    gbr_model = gbr_models[country]
    y_pred_gbr = gbr_model.predict(X_test)
    r_squared_gbr = r2_score(y_test, y_pred_gbr)
    
    # Append the summary to the list
    gbr_model_summary.append({
        'Country': country,
        'Gradient Boosting R^2 Score': r_squared_gbr
    })

# Create a DataFrame to display the Gradient Boosting model accuracy summary
gbr_model_summary_df = pd.DataFrame(gbr_model_summary)

gbr_model_summary_df

# Prepare a list to store the summary of the Gradient Boosting model's 
# intercept (using the mean), coefficient approximation (feature importance), and R^2 for each country

gbr_summary_list = []  # List to store summaries for each country

for country, gbr_model in gbr_models.items():
    # Calculate the mean of predictions as an approximation for the intercept
    mean_intercept = np.mean(gbr_model.predict(X_train))  # Assuming X_train contains the historical features

    # Get feature importances from the Gradient Boosting model
    feature_importance = gbr_model.feature_importances_

    # Calculate the R^2 score of the model
    r_squared = gbr_model.score(X_train, y_train)

    # Prepare a dictionary for the current country's summary
    country_summary = {
        'Country': country,
        'Mean_Intercept': mean_intercept,
        'Feature_Importance': feature_importance,
        'R_squared': r_squared
    }

    # Append the dictionary to the list
    gbr_summary_list.append(country_summary)

gbr_summary_df = pd.DataFrame(gbr_summary_list)

gbr_summary_df

prophet_data = data[['Country', 'Year', 'Expenditure_Constant_2022']].copy() 
USA_df = prophet_data[prophet_data['Country'] == 'United States of America']

#set up country dataframe to be used by prophet
USA_df2= USA_df[['Year', 'Expenditure_Constant_2022']]

# Rename the columns to 'ds' and 'y' for use by Prophet and convert time date time stamp expected by Prophet
USA_df2 = USA_df2.rename(columns={'Year': 'ds', 'Expenditure_Constant_2022': 'y'})
USA_df2['ds'] = to_datetime(USA_df2['ds'], format='%Y')

# Fit the model by instantiating a new Prophet object
m = Prophet()

#Call the fit method and pass in the historical dataframe. Fitting should take 1-5 seconds.m.fit(df)
m.fit(USA_df2)
future = m.make_future_dataframe(periods=10, freq = "Y", include_history= True)
forecast = m.predict(future)
fig1 = m.plot(forecast)
plt.title('USA Military Expenditure Forecast')
plt.xlabel('Year')
plt.ylabel('Expenditure_Constant_2022')
plt.show()

21:21:31 - cmdstanpy - INFO - Chain [1] start processing
21:21:31 - cmdstanpy - INFO - Chain [1] done processing
C:\Users\swapn\AppData\Roaming\Python\Python312\site-packages\prophet\forecaster.py:1854: FutureWarning: 'Y' is deprecated and will be removed in a future version, please use 'YE' instead.
  dates = pd.date_range(

# Perform cross-validation
df_cv = cross_validation(m, initial='730 days', period='180 days', horizon='365 days')

# Compute performance metrics
df_p = performance_metrics(df_cv)
print(df_p)

  0%|          | 0/41 [00:00<?, ?it/s]

21:21:32 - cmdstanpy - INFO - Chain [1] start processing
21:21:32 - cmdstanpy - INFO - Chain [1] done processing
21:21:33 - cmdstanpy - INFO - Chain [1] start processing
21:21:33 - cmdstanpy - INFO - Chain [1] done processing
21:21:34 - cmdstanpy - INFO - Chain [1] start processing
21:21:34 - cmdstanpy - INFO - Chain [1] done processing
21:21:34 - cmdstanpy - INFO - Chain [1] start processing
21:21:35 - cmdstanpy - INFO - Chain [1] done processing
21:21:35 - cmdstanpy - INFO - Chain [1] start processing
21:21:45 - cmdstanpy - INFO - Chain [1] done processing
21:21:45 - cmdstanpy - INFO - Chain [1] start processing
21:21:55 - cmdstanpy - INFO - Chain [1] done processing
21:21:55 - cmdstanpy - INFO - Chain [1] start processing
21:22:11 - cmdstanpy - INFO - Chain [1] done processing
21:22:12 - cmdstanpy - INFO - Chain [1] start processing
21:22:28 - cmdstanpy - INFO - Chain [1] done processing
21:22:28 - cmdstanpy - INFO - Chain [1] start processing
21:22:28 - cmdstanpy - INFO - Chain [1] done processing
21:22:28 - cmdstanpy - INFO - Chain [1] start processing
21:22:29 - cmdstanpy - INFO - Chain [1] done processing
21:22:29 - cmdstanpy - INFO - Chain [1] start processing
21:22:29 - cmdstanpy - INFO - Chain [1] done processing
21:22:29 - cmdstanpy - INFO - Chain [1] start processing
21:22:29 - cmdstanpy - INFO - Chain [1] done processing
21:22:29 - cmdstanpy - INFO - Chain [1] start processing
21:22:29 - cmdstanpy - INFO - Chain [1] done processing
21:22:30 - cmdstanpy - INFO - Chain [1] start processing
21:22:30 - cmdstanpy - INFO - Chain [1] done processing
21:22:30 - cmdstanpy - INFO - Chain [1] start processing
21:22:30 - cmdstanpy - INFO - Chain [1] done processing
21:22:30 - cmdstanpy - INFO - Chain [1] start processing
21:22:30 - cmdstanpy - INFO - Chain [1] done processing
21:22:31 - cmdstanpy - INFO - Chain [1] start processing
21:22:31 - cmdstanpy - INFO - Chain [1] done processing
21:22:31 - cmdstanpy - INFO - Chain [1] start processing
21:22:31 - cmdstanpy - INFO - Chain [1] done processing
21:22:31 - cmdstanpy - INFO - Chain [1] start processing
21:22:32 - cmdstanpy - INFO - Chain [1] done processing
21:22:32 - cmdstanpy - INFO - Chain [1] start processing
21:22:32 - cmdstanpy - INFO - Chain [1] done processing
21:22:32 - cmdstanpy - INFO - Chain [1] start processing
21:22:33 - cmdstanpy - INFO - Chain [1] done processing
21:22:33 - cmdstanpy - INFO - Chain [1] start processing
21:22:33 - cmdstanpy - INFO - Chain [1] done processing
21:22:33 - cmdstanpy - INFO - Chain [1] start processing
21:22:33 - cmdstanpy - INFO - Chain [1] done processing
21:22:34 - cmdstanpy - INFO - Chain [1] start processing
21:22:34 - cmdstanpy - INFO - Chain [1] done processing
21:22:34 - cmdstanpy - INFO - Chain [1] start processing
21:22:34 - cmdstanpy - INFO - Chain [1] done processing
21:22:34 - cmdstanpy - INFO - Chain [1] start processing
21:22:35 - cmdstanpy - INFO - Chain [1] done processing
21:22:35 - cmdstanpy - INFO - Chain [1] start processing
21:22:35 - cmdstanpy - INFO - Chain [1] done processing
21:22:35 - cmdstanpy - INFO - Chain [1] start processing
21:22:35 - cmdstanpy - INFO - Chain [1] done processing
21:22:35 - cmdstanpy - INFO - Chain [1] start processing
21:22:36 - cmdstanpy - INFO - Chain [1] done processing
21:22:36 - cmdstanpy - INFO - Chain [1] start processing
21:22:36 - cmdstanpy - INFO - Chain [1] done processing
21:22:36 - cmdstanpy - INFO - Chain [1] start processing
21:22:37 - cmdstanpy - INFO - Chain [1] done processing
21:22:37 - cmdstanpy - INFO - Chain [1] start processing
21:22:37 - cmdstanpy - INFO - Chain [1] done processing
21:22:37 - cmdstanpy - INFO - Chain [1] start processing
21:22:38 - cmdstanpy - INFO - Chain [1] done processing
21:22:38 - cmdstanpy - INFO - Chain [1] start processing
21:22:38 - cmdstanpy - INFO - Chain [1] done processing
21:22:38 - cmdstanpy - INFO - Chain [1] start processing
21:22:38 - cmdstanpy - INFO - Chain [1] done processing
21:22:39 - cmdstanpy - INFO - Chain [1] start processing
21:22:39 - cmdstanpy - INFO - Chain [1] done processing
21:22:39 - cmdstanpy - INFO - Chain [1] start processing
21:22:39 - cmdstanpy - INFO - Chain [1] done processing
21:22:40 - cmdstanpy - INFO - Chain [1] start processing
21:22:40 - cmdstanpy - INFO - Chain [1] done processing
21:22:40 - cmdstanpy - INFO - Chain [1] start processing
21:22:40 - cmdstanpy - INFO - Chain [1] done processing
21:22:40 - cmdstanpy - INFO - Chain [1] start processing
21:22:41 - cmdstanpy - INFO - Chain [1] done processing
21:22:41 - cmdstanpy - INFO - Chain [1] start processing
21:22:41 - cmdstanpy - INFO - Chain [1] done processing

    horizon           mse           rmse            mae      mape     mdape  \
0   96 days  2.502831e+09   50028.298004   42244.965449  0.053806  0.042097   
1  101 days  3.682110e+09   60680.393849   55675.224043  0.068985  0.068527   
2  106 days  3.790948e+09   61570.675560   56961.494120  0.068548  0.067655   
3  112 days  3.728111e+09   61058.261477   55593.885857  0.065759  0.067655   
4  117 days  1.857993e+09   43104.447679   34160.981158  0.039274  0.035645   
5  122 days  1.479856e+09   38468.893014   31249.614687  0.033115  0.035645   
6  127 days  6.822388e+09   82597.748197   58259.447889  0.061904  0.040709   
7  133 days  1.981269e+10  140757.544082  110411.593548  0.123698  0.112915   
8  138 days  3.495308e+10  186957.418691  171744.671410  0.200333  0.216810   
9  143 days  4.817319e+10  219483.916647  216323.762017  0.261079  0.285790   
10 148 days  5.305978e+10  230347.086701  229893.225657  0.286113  0.285790   
11 154 days  4.961047e+10  222734.076939  221781.933117  0.282652  0.284475   
12 159 days  3.922574e+10  198054.888516  194739.664486  0.249175  0.259956   
13 164 days  2.677118e+10  163619.022890  155625.109374  0.197411  0.213874   
14 169 days  1.637141e+10  127950.814051  112948.629230  0.141384  0.134865   
15 175 days  7.229441e+09   85026.118691   75044.377718  0.090655  0.073881   
16 180 days  3.374434e+09   58089.876664   55569.852989  0.064741  0.060587   
17 260 days  1.807733e+09   42517.441180   39586.425370  0.046738  0.046445   
18 265 days  1.869546e+09   43238.244484   40397.338595  0.049378  0.051726   
19 271 days  1.533366e+09   39158.213062   35628.345252  0.044564  0.042097   
20 276 days  2.502831e+09   50028.298004   42244.965449  0.053806  0.042097   
21 281 days  3.682110e+09   60680.393849   55675.224043  0.068985  0.068527   
22 286 days  3.790948e+09   61570.675560   56961.494120  0.068548  0.067655   
23 292 days  3.728111e+09   61058.261477   55593.885857  0.065759  0.067655   
24 297 days  1.857993e+09   43104.447679   34160.981158  0.039274  0.035645   
25 302 days  1.479856e+09   38468.893014   31249.614687  0.033115  0.035645   
26 307 days  6.822388e+09   82597.748197   58259.447889  0.061904  0.040709   
27 313 days  1.981269e+10  140757.544082  110411.593548  0.123698  0.112915   
28 318 days  3.495308e+10  186957.418691  171744.671410  0.200333  0.216810   
29 323 days  4.817319e+10  219483.916647  216323.762017  0.261079  0.285790   
30 328 days  5.305978e+10  230347.086701  229893.225657  0.286113  0.285790   
31 334 days  4.961047e+10  222734.076939  221781.933117  0.282652  0.284475   
32 339 days  3.922574e+10  198054.888516  194739.664486  0.249175  0.259956   
33 344 days  2.677118e+10  163619.022890  155625.109374  0.197411  0.213874   
34 349 days  1.637141e+10  127950.814051  112948.629230  0.141384  0.134865   
35 355 days  7.229441e+09   85026.118691   75044.377718  0.090655  0.073881   
36 360 days  3.374434e+09   58089.876664   55569.852989  0.064741  0.060587   
37 365 days  2.198867e+09   46892.075782   46093.054930  0.052893  0.050186   

       smape  coverage  
0   0.052331      0.00  
1   0.066568      0.00  
2   0.066153      0.00  
3   0.063288      0.25  
4   0.038154      0.50  
5   0.032420      0.50  
6   0.058357      0.50  
7   0.112168      0.25  
8   0.178570      0.00  
9   0.229769      0.00  
10  0.250181      0.00  
11  0.247473      0.00  
12  0.220770      0.00  
13  0.177889      0.25  
14  0.129434      0.50  
15  0.085557      0.75  
16  0.062514      1.00  
17  0.045692      0.75  
18  0.048213      0.50  
19  0.043853      0.25  
20  0.052331      0.00  
21  0.066568      0.00  
22  0.066153      0.00  
23  0.063288      0.25  
24  0.038154      0.50  
25  0.032420      0.50  
26  0.058357      0.50  
27  0.112168      0.25  
28  0.178570      0.00  
29  0.229769      0.00  
30  0.250181      0.00  
31  0.247473      0.00  
32  0.220770      0.00  
33  0.177889      0.25  
34  0.129434      0.50  
35  0.085557      0.75  
36  0.062514      1.00  
37  0.051480      1.00

# Selecting relevant features for clustering
features_for_clustering = [
    "GovtSpending", 
    "GDP", 
    "PerCapita", 
    "ExpenditureConstant",
    "TotalScore",
    "Deadliness", 
    #"Diffusion", 
    #"Danger", 
    #"Fragmentation"
]

# Extracting the data
clustering_data = data_for_clustering[features_for_clustering]

# Normalizing the data using StandardScaler
scaler = StandardScaler()
normalized_data = scaler.fit_transform(clustering_data)

# Converting the normalized data back to a DataFrame for easier inspection
normalized_df = pd.DataFrame(normalized_data, columns=features_for_clustering)

# Displaying the first few rows of the normalized data
normalized_df.head()

# Applying PCA to reduce dimensions
pca = PCA(n_components=3)  # Reduce to 3 dimensions for faster computation and visualization
reduced_data = pca.fit_transform(normalized_data)

# Check the explained variance ratio to ensure sufficient information is retained
explained_variance_ratio = pca.explained_variance_ratio_
explained_variance_ratio.sum()

print("Explained Variance Ratio:", explained_variance_ratio)

# PCA components (eigenvectors)
components = pca.components_

# Create a DataFrame for feature contributions
feature_contribution = pd.DataFrame(components, columns=features_for_clustering, 
                                     index=[f'PC{i+1}' for i in range(components.shape[0])])

print("Feature Contributions to PCA Components:")
print(feature_contribution)

Explained Variance Ratio: [0.35639507 0.3023877  0.18063604]
Feature Contributions to PCA Components:
     GovtSpending       GDP  PerCapita  ExpenditureConstant  TotalScore  \
PC1      0.532087  0.558538   0.352505             0.209272    0.336779   
PC2     -0.151693 -0.213752  -0.411999            -0.275108    0.589410   
PC3     -0.418062 -0.314193   0.337101             0.739804    0.190152   

     Deadliness  
PC1    0.351346  
PC2    0.581780  
PC3    0.171469

# Apply Kernel PCA with a polynomial kernel
n_components = 2  # Number of components to retain
kernel_pca = KernelPCA(n_components=n_components, kernel='poly', degree=3, gamma=1, coef0=1)  
# Degree and coef0 can be tuned
reduced_data = kernel_pca.fit_transform(normalized_data)

# Perform KMeans clustering on the Kernel PCA-reduced data
kmeans = KMeans(n_clusters=2, random_state=42)
poly_labels = kmeans.fit_predict(reduced_data)

# Calculate Silhouette Score
poly_silhouette = silhouette_score(reduced_data, poly_labels)
print(f"Silhouette Score with Polynomial Kernel PCA: {poly_silhouette:.4f}")

# Visualize the clusters from Polynomial Kernel PCA
plt.figure(figsize=(10, 6))
plt.scatter(reduced_data[:, 0], reduced_data[:, 1], c=poly_labels, cmap='viridis', s=50, alpha=0.7)
plt.title(f'KMeans Clustering with Polynomial Kernel PCA ({n_components} Clusters)')
plt.xlabel('Kernel Principal Component 1')
plt.ylabel('Kernel Principal Component 2')
plt.colorbar(label='Cluster')
plt.grid(True)
plt.show()

Silhouette Score with Polynomial Kernel PCA: 0.9887

# Elbow method to determine the optimal number of clusters using reduced data
inertia_reduced = []
range_n_clusters = range(1, 11)

for k in range_n_clusters:
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(reduced_data)
    inertia_reduced.append(kmeans.inertia_)

# Find the elbow point
knee_locator = KneeLocator(range_n_clusters, inertia_reduced, curve="convex", direction="decreasing")
optimal_k = knee_locator.knee

# Plotting the Elbow curve
plt.figure(figsize=(10, 6))
plt.plot(range_n_clusters, inertia_reduced, marker='o', linestyle='--', color='blue', label='Inertia')
plt.title('Elbow Method Analysis (PCA-Reduced Data)', fontsize=14, fontweight='bold')
plt.xlabel('Number of Clusters (k)', fontsize=12)
plt.ylabel('Inertia', fontsize=12)
plt.xticks(range_n_clusters, fontsize=10)
plt.yticks(fontsize=10)

# Highlight the elbow point
plt.scatter(optimal_k, inertia_reduced[optimal_k - 1], color='red', s=100, zorder=5, label=f'Optimal k = {optimal_k}')
plt.axvline(x=optimal_k, linestyle='--', color='red')
plt.annotate(f'Optimal k = {optimal_k}', 
             xy=(optimal_k, inertia_reduced[optimal_k - 1]), 
             xytext=(optimal_k + 1, inertia_reduced[optimal_k - 1] + 1000), 
             arrowprops=dict(facecolor='black', arrowstyle='->'),
             fontsize=10)

plt.grid(color='gray', linestyle='--', linewidth=0.5, alpha=0.7)
plt.legend(fontsize=10)
plt.tight_layout()
plt.show()

# selecting the optimal number of clusters
optimal_clusters = 4

def perform_kmeans_clustering(data, n_clusters, print_metrics=True, visual=True):

    # Perform KMeans clustering with the chosen number of clusters
    kmeans = KMeans(n_clusters=n_clusters, random_state=42)
    labels = kmeans.fit_predict(data)

    # Calculate clustering evaluation metrics
    silhouette = silhouette_score(data, labels)
    calinski_harabasz = calinski_harabasz_score(data, labels)
    davies_bouldin = davies_bouldin_score(data, labels)

    # Store the evaluation metrics in a dictionary for easier reporting
    metrics = {
        "Silhouette Score": silhouette,
        "Calinski-Harabasz Score": calinski_harabasz,
        "Davies-Bouldin Score": davies_bouldin
    }

    # Print the K-Means evaluation metrics
    if print_metrics:
        print("KMeans Clustering Metrics")
        for metric, value in metrics.items():
            print(f"{metric}: {value:.4f}")

    # Visualize the clusters
    if visual:
        plt.figure(figsize=(10, 6))
        plt.scatter(data[:, 0], data[:, 1], c=labels, cmap='viridis', s=50, alpha=0.7)
        plt.title('KMeans Clustering (PCA-Reduced Data)')
        plt.xlabel('Principal Component 1')
        plt.ylabel('Principal Component 2')
        plt.colorbar(label='Cluster')
        plt.grid(True)
        plt.show()

    return metrics, labels

# Example usage
kmeans_metrics, final_labels = perform_kmeans_clustering(reduced_data, optimal_clusters, print_metrics=True, visual=True)

KMeans Clustering Metrics
Silhouette Score: 0.9919
Calinski-Harabasz Score: 24071.5269
Davies-Bouldin Score: 0.1980

def perform_gmm_clustering(data, n_components,  print_metrics=True, visual=True):
    # Run Gaussian Mixture Model (GMM) clustering
    gmm = GaussianMixture(n_components=n_components, random_state=42)
    gmm_labels = gmm.fit_predict(data)

    # Calculate evaluation metrics for GMM clustering
    silhouette = silhouette_score(data, gmm_labels)
    calinski_harabasz = calinski_harabasz_score(data, gmm_labels)
    davies_bouldin = davies_bouldin_score(data, gmm_labels)

    # Store the evaluation metrics in a dictionary
    gmm_metrics = {
        "Silhouette Score": silhouette,
        "Calinski-Harabasz Score": calinski_harabasz,
        "Davies-Bouldin Score": davies_bouldin
    }

    # Print the evaluation metrics for GMM clustering
    if print_metrics:
        print("GMM Evaluation Metrics")
        for metric, value in gmm_metrics.items():
            print(f"{metric}: {value:.4f}")

    # Visualize GMM clusters
    if visual:
        plt.figure(figsize=(10, 6))
        plt.scatter(data[:, 0], data[:, 1], c=gmm_labels, cmap='viridis', s=50, alpha=0.7)
        plt.title('Gaussian Mixture Model Clustering (PCA-Reduced Data)')
        plt.xlabel('Principal Component 1')
        plt.ylabel('Principal Component 2')
        plt.colorbar(label='Cluster')
        plt.grid(True)
        plt.show()

    return gmm_metrics, gmm_labels

# Example usage
gmm_metrics, gmm_labels = perform_gmm_clustering(reduced_data, optimal_clusters,  print_metrics=True, visual=True)

GMM Evaluation Metrics
Silhouette Score: 0.9821
Calinski-Harabasz Score: 7458.1424
Davies-Bouldin Score: 0.4536

def perform_agglomerative_clustering(data, n_clusters, print_metrics=True, visual=True):
    # Perform Agglomerative Clustering
    agg_clustering = AgglomerativeClustering(n_clusters=n_clusters, linkage='ward')
    agg_labels = agg_clustering.fit_predict(data)

    # Compute linkage matrix for dendrogram visualization
    linkage_matrix = linkage(data, method='ward')

    # Compute clustering evaluation metrics
    silhouette = silhouette_score(data, agg_labels)
    calinski_harabasz = calinski_harabasz_score(data, agg_labels)
    davies_bouldin = davies_bouldin_score(data, agg_labels)

    # Store the evaluation metrics in a dictionary
    agg_metrics = {
        "Silhouette Score": silhouette,
        "Calinski-Harabasz Score": calinski_harabasz,
        "Davies-Bouldin Score": davies_bouldin
    }

    # Print the evaluation metrics for Agglomerative Clustering
    if print_metrics:
        print("Agglomerative Clustering Evaluation Metrics")
        for metric, value in agg_metrics.items():
            print(f"{metric}: {value:.4f}")

    # Plot dendrogram for hierarchical clustering visualization
    if visual:
        plt.figure(figsize=(12, 8))
        dendrogram(linkage_matrix, truncate_mode='level', p=5, color_threshold=0.5 * max(linkage_matrix[:, 2]))
        plt.title("Hierarchical Clustering Dendrogram")
        plt.xlabel("Sample Index or Cluster Size")
        plt.ylabel("Distance")
        plt.grid(True)
        plt.show()

    return agg_metrics, agg_labels

# Example usage
agg_metrics, agg_labels = perform_agglomerative_clustering(reduced_data, optimal_clusters, print_metrics=True, visual=True)

Agglomerative Clustering Evaluation Metrics
Silhouette Score: 0.9919
Calinski-Harabasz Score: 23860.9024
Davies-Bouldin Score: 0.1916

def perform_dbscan_clustering(data, eps, min_samples=5, print_metrics=True, visual=True):
    # Run DBSCAN clustering
    dbscan = DBSCAN(eps=eps, min_samples=min_samples)
    dbscan_labels = dbscan.fit_predict(data)

    # Calculate evaluation metrics for DBSCAN (only for clusters with more than one label)
    if len(set(dbscan_labels)) > 1:
        silhouette = silhouette_score(data, dbscan_labels)
        calinski_harabasz = calinski_harabasz_score(data, dbscan_labels)
        davies_bouldin = davies_bouldin_score(data, dbscan_labels)
    else:
        silhouette, calinski_harabasz, davies_bouldin = None, None, None

    # Store the evaluation metrics in a dictionary
    dbscan_metrics = {
        "Silhouette Score": silhouette,
        "Calinski-Harabasz Score": calinski_harabasz,
        "Davies-Bouldin Score": davies_bouldin
    }

    # Print the evaluation metrics
    if print_metrics:
        print("DBSCAN Evaluation Metrics: Epsilon =", eps)
        for metric, value in dbscan_metrics.items():
            print(f"{metric}: {value:.4f}" if value is not None else f"{metric}: No clusters found")

    # Visualize DBSCAN clusters
    if visual:
        plt.figure(figsize=(10, 6))
        plt.scatter(data[:, 0], data[:, 1], c=dbscan_labels, cmap='viridis', s=50, alpha=0.7)
        plt.title('DBSCAN Clustering')
        plt.xlabel('Principal Component 1')
        plt.ylabel('Principal Component 2')
        plt.colorbar(label='Cluster')
        plt.grid(True)
        plt.show()

    return dbscan_metrics, dbscan_labels

optimal_eps = 0.5 # optimal value for DBSCAN's epsilon parameter
dbscan_metrics, dbscan_labels = perform_dbscan_clustering(reduced_data, optimal_eps, print_metrics=True, visual=True)

DBSCAN Evaluation Metrics: Epsilon = 0.5
Silhouette Score: 0.5104
Calinski-Harabasz Score: 50.0124
Davies-Bouldin Score: 1.9200

# Define the range of clusters and eps values to test
cluster_range = [2, 3, 4, 5, 6, 7, 8, 9]
eps_range = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]

# Initialize dictionaries to store metrics for each clustering method
kmeans_metrics_dict = {}
gmm_metrics_dict = {}
agg_metrics_dict = {}
dbscan_metrics_dict = {}

# Run clustering for each cluster count and store the metrics
for n_clusters in cluster_range:
    kmeans_metrics, _ = perform_kmeans_clustering(reduced_data, n_clusters, print_metrics=False, visual=False)
    kmeans_metrics_dict[n_clusters] = kmeans_metrics
    gmm_metrics, _ = perform_gmm_clustering(reduced_data, n_clusters, print_metrics=False, visual=False)
    gmm_metrics_dict[n_clusters] = gmm_metrics
    agg_metrics, _ = perform_agglomerative_clustering(reduced_data, n_clusters, print_metrics=False, visual=False)
    agg_metrics_dict[n_clusters] = agg_metrics

# Run DBSCAN clustering for each eps value and store the metrics
for eps in eps_range:
    dbscan_metrics, _ = perform_dbscan_clustering(reduced_data, eps, print_metrics=False, visual=False)
    dbscan_metrics_dict[eps] = dbscan_metrics

# Convert the metrics dictionaries to DataFrames for easier visualization
kmeans_metrics_df = pd.DataFrame(kmeans_metrics_dict).T
gmm_metrics_df = pd.DataFrame(gmm_metrics_dict).T
agg_metrics_df = pd.DataFrame(agg_metrics_dict).T
dbscan_metrics_df = pd.DataFrame(dbscan_metrics_dict).T

# Plot the metrics for each clustering method
fig, axes = plt.subplots(4, 3, figsize=(18, 20), sharex=True)

# Plot metrics
def plot_metrics(axes, cluster_range, metrics_df, algo_name, score_name):
    axes.plot(cluster_range, metrics_df[score_name], marker='o', linestyle='--', color='b')
    axes.set_title(f'{algo_name}-{score_name}')  
    axes.set_ylabel('Score')
    if algo_name == 'DBSCAN':
        axes.set_xlabel('Epsilon')
    else:
        axes.set_xlabel('Clusters')

plot_metrics(axes[0,0], cluster_range, kmeans_metrics_df, 'KMeans', 'Silhouette Score')
plot_metrics(axes[0,1], cluster_range, kmeans_metrics_df, 'KMeans', 'Calinski-Harabasz Score')
plot_metrics(axes[0,2], cluster_range, kmeans_metrics_df, 'KMeans', 'Davies-Bouldin Score')
plot_metrics(axes[1,0], cluster_range, gmm_metrics_df, 'GMM', 'Silhouette Score')
plot_metrics(axes[1,1], cluster_range, gmm_metrics_df, 'GMM', 'Calinski-Harabasz Score')
plot_metrics(axes[1,2], cluster_range, gmm_metrics_df, 'GMM', 'Davies-Bouldin Score')
plot_metrics(axes[2,0], cluster_range, agg_metrics_df, 'Agglomerative', 'Silhouette Score')
plot_metrics(axes[2,1], cluster_range, agg_metrics_df, 'Agglomerative', 'Calinski-Harabasz Score')
plot_metrics(axes[2,2], cluster_range, agg_metrics_df, 'Agglomerative', 'Davies-Bouldin Score')
plot_metrics(axes[3,0], eps_range, dbscan_metrics_df, 'DBSCAN', 'Silhouette Score')
plot_metrics(axes[3,1], eps_range, dbscan_metrics_df, 'DBSCAN', 'Calinski-Harabasz Score')
plot_metrics(axes[3,2], eps_range, dbscan_metrics_df, 'DBSCAN', 'Davies-Bouldin Score')
   
plt.tight_layout()
plt.show()

# print the metrics for each clustering method in a tabular format
print("KMeans Metrics:")
print(kmeans_metrics_df)
print("\nGMM Metrics:")
print(gmm_metrics_df)
print("\nAgglomerative Metrics:")
print(agg_metrics_df)
print("\nDBSCAN Metrics:")
print(dbscan_metrics_df)

KMeans Metrics:
   Silhouette Score  Calinski-Harabasz Score  Davies-Bouldin Score
2          0.988711              1495.739776              0.240315
3          0.992138              3668.948110              0.136262
4          0.991880             24071.526904              0.197978
5          0.990031             26941.894976              0.301172
6          0.989663             23976.677108              0.224383
7          0.987734             36354.252282              0.214310
8          0.987958             41025.957959              0.179783
9          0.961958             62394.518810              0.232444

GMM Metrics:
   Silhouette Score  Calinski-Harabasz Score  Davies-Bouldin Score
2          0.977369               892.234837              0.702360
3          0.992138              3668.948110              0.136262
4          0.982099              7458.142415              0.453593
5          0.977323             11589.768293              0.559008
6          0.976955              9731.567537              0.462436
7          0.978502             10123.726658              0.479276
8          0.977841              9395.064499              0.459485
9          0.604006              9126.578867              0.603746

Agglomerative Metrics:
   Silhouette Score  Calinski-Harabasz Score  Davies-Bouldin Score
2          0.996496              4138.303415              0.002559
3          0.990418              8109.212032              0.159652
4          0.991901             23860.902391              0.191648
5          0.989906             24946.192382              0.298318
6          0.987638             32150.406775              0.295904
7          0.983638             37075.601426              0.325685
8          0.952023             44911.779048              0.391883
9          0.952253             55927.215660              0.352276

DBSCAN Metrics:
     Silhouette Score  Calinski-Harabasz Score  Davies-Bouldin Score
0.2          0.457582                13.109143              1.965495
0.3          0.546640                21.936386              1.958826
0.4          0.534797                34.685951              1.941840
0.5          0.510430                50.012376              1.920019
0.6          0.762987                53.404171              1.929049
0.7          0.793035                80.108842              1.929897
0.8          0.643066                34.248305              1.937013
0.9          0.647645                29.987315              1.926479

# Remove non-numeric columns for the feature set preparation
numeric_columns = data_cleaned.select_dtypes(include=[np.number]).columns

# Feature selection and data preparation
features = [
    'Year', 'index_level', 'total_rank_from_avg_rank', 'total_score_rank',
    'Deadliness_raw', 'Diffusion_raw', 'Danger_raw', 'Fragmentation_raw',
    'Expenditure-Share_of_Govt_spending', 'Expenditure-Share_of_GDP', 'Expenditure_Per_Capita'
]
features = [feature for feature in features if feature in numeric_columns]

target = 'Expenditure_Constant_2022'

# Prepare the feature set and target vector
data_for_model = data_cleaned[['Country', 'Year', target] + features].dropna()

# Convert the 'Year' column to integers to avoid any data type inconsistencies
data_for_model['Year'] = data_for_model['Year'].astype(int)

# Remove any rows with missing 'Year' values to ensure consistency
data_for_model = data_for_model.dropna(subset=['Year'])

# Extract the latest year for each country using groupby and then merge to get the latest expenditure data
latest_year_per_country = data_for_model.groupby('Country')['Year'].max().reset_index()

# Drop extra 'Year' columns to avoid duplicates
data_for_merge = data_for_model.drop(columns=['Year'])

# Merge back to get data for the latest year for each country
latest_expenditure = pd.merge(data_for_merge, latest_year_per_country, on=['Country'], how='inner')

# Get the top 15 countries based on the latest expenditure values
top_15_countries = latest_expenditure.groupby('Country')[target].max().nlargest(15).index
# Filter the data for the top 15 countries
top_15_data = data_for_model[data_for_model['Country'].isin(top_15_countries)]

# Prepare feature matrix (X) and target vector (y)
X = top_15_data[features]
y = top_15_data[target]

# Standardize features
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Split the data into training and testing sets (80% train, 20% test)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)

# Defining the parameter grid for Random Forest
param_grid = {
    'n_estimators': [100, 200, 300],
    'max_depth': [10, 20, None],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4],
    'max_features': ['auto', 'sqrt', 'log2']
}

# Initialize the Random Forest Regressor
rf_model = RandomForestRegressor(random_state=42)

# Setting up the GridSearchCV
grid_search = GridSearchCV(estimator=rf_model,
                           param_grid=param_grid,
                           cv=3,
                           n_jobs=-1,
                           verbose=2,
                           scoring='r2')

# Fitting GridSearchCV on the training data
grid_search.fit(X_train, y_train)

# Displaying the best parameters from the grid search
best_params = grid_search.best_params_
print("Best parameters found: ", best_params)

# Using the best estimator to evaluate on the test data
best_rf_model = grid_search.best_estimator_
r2_score = best_rf_model.score(X_test, y_test)
print(f"R^2 Score for the best Random Forest model: {r2_score}")

Fitting 3 folds for each of 243 candidates, totalling 729 fits

c:\ProgramData\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:547: FitFailedWarning: 
243 fits failed out of a total of 729.
The score on these train-test partitions for these parameters will be set to nan.
If these failures are not expected, you can try to debug them by setting error_score='raise'.

Below are more details about the failures:
--------------------------------------------------------------------------------
185 fits failed with the following error:
Traceback (most recent call last):
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py", line 895, in _fit_and_score
    estimator.fit(X_train, y_train, **fit_params)
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\base.py", line 1467, in wrapper
    estimator._validate_params()
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\base.py", line 666, in _validate_params
    validate_parameter_constraints(
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\utils\_param_validation.py", line 95, in validate_parameter_constraints
    raise InvalidParameterError(
sklearn.utils._param_validation.InvalidParameterError: The 'max_features' parameter of RandomForestRegressor must be an int in the range [1, inf), a float in the range (0.0, 1.0], a str among {'log2', 'sqrt'} or None. Got 'auto' instead.

--------------------------------------------------------------------------------
58 fits failed with the following error:
Traceback (most recent call last):
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py", line 895, in _fit_and_score
    estimator.fit(X_train, y_train, **fit_params)
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\base.py", line 1467, in wrapper
    estimator._validate_params()
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\base.py", line 666, in _validate_params
    validate_parameter_constraints(
  File "c:\ProgramData\anaconda3\Lib\site-packages\sklearn\utils\_param_validation.py", line 95, in validate_parameter_constraints
    raise InvalidParameterError(
sklearn.utils._param_validation.InvalidParameterError: The 'max_features' parameter of RandomForestRegressor must be an int in the range [1, inf), a float in the range (0.0, 1.0], a str among {'sqrt', 'log2'} or None. Got 'auto' instead.

  warnings.warn(some_fits_failed_message, FitFailedWarning)
c:\ProgramData\anaconda3\Lib\site-packages\sklearn\model_selection\_search.py:1051: UserWarning: One or more of the test scores are non-finite: [       nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan 0.99147593 0.99067179 0.99037055
 0.98980537 0.98684391 0.98710208 0.96519808 0.96089731 0.96150901
 0.98373235 0.97995591 0.98140769 0.98342261 0.9796001  0.97978418
 0.96257031 0.95859151 0.96021209 0.96315601 0.95461746 0.95682406
 0.96315601 0.95461746 0.95682406 0.9580552  0.95133184 0.95346238
 0.99147593 0.99067179 0.99037055 0.98980537 0.98684391 0.98710208
 0.96519808 0.96089731 0.96150901 0.98373235 0.97995591 0.98140769
 0.98342261 0.9796001  0.97978418 0.96257031 0.95859151 0.96021209
 0.96315601 0.95461746 0.95682406 0.96315601 0.95461746 0.95682406
 0.9580552  0.95133184 0.95346238        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
 0.98993773 0.98791543 0.98898602 0.98891951 0.98564412 0.98654724
 0.96469892 0.96042635 0.96108508 0.98463325 0.9802256  0.98081682
 0.98329412 0.97952853 0.97952191 0.96250046 0.95857727 0.96025448
 0.96306094 0.95445917 0.95672779 0.96306094 0.95445917 0.95672779
 0.9580552  0.95133184 0.95336487 0.98993773 0.98791543 0.98898602
 0.98891951 0.98564412 0.98654724 0.96469892 0.96042635 0.96108508
 0.98463325 0.9802256  0.98081682 0.98329412 0.97952853 0.97952191
 0.96250046 0.95857727 0.96025448 0.96306094 0.95445917 0.95672779
 0.96306094 0.95445917 0.95672779 0.9580552  0.95133184 0.95336487
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan        nan        nan        nan
        nan        nan        nan 0.98993773 0.98791543 0.98898602
 0.98891951 0.98564412 0.98654724 0.96469892 0.96042635 0.96108508
 0.98463325 0.9802256  0.98081682 0.98329412 0.97952853 0.97952191
 0.96250046 0.95857727 0.96025448 0.96306094 0.95445917 0.95672779
 0.96306094 0.95445917 0.95672779 0.9580552  0.95133184 0.95336487
 0.98993773 0.98791543 0.98898602 0.98891951 0.98564412 0.98654724
 0.96469892 0.96042635 0.96108508 0.98463325 0.9802256  0.98081682
 0.98329412 0.97952853 0.97952191 0.96250046 0.95857727 0.96025448
 0.96306094 0.95445917 0.95672779 0.96306094 0.95445917 0.95672779
 0.9580552  0.95133184 0.95336487]
  warnings.warn(

Best parameters found:  {'max_depth': 10, 'max_features': 'sqrt', 'min_samples_leaf': 1, 'min_samples_split': 2, 'n_estimators': 100}
R^2 Score for the best Random Forest model: 0.9955922371327366

# Initialize the Random Forest Regressor with the tuned parameters
rf_model = RandomForestRegressor(
    n_estimators=100, 
    max_depth=10, 
    max_features='sqrt', 
    min_samples_leaf=1, 
    min_samples_split=2, 
    random_state=42
)

# Fit the Random Forest model on the training data
rf_model.fit(X_train, y_train)

# Evaluate the model on the test data
r2_score = rf_model.score(X_test, y_test)
print(f"R^2 Score for the tuned Random Forest model: {r2_score}")

R^2 Score for the tuned Random Forest model: 0.9955922371327366

# Feature Importance Analysis (for Random Forest)
feature_importances = pd.Series(rf_model.feature_importances_, index=X.columns)
feature_importances.nlargest(10).plot(kind='barh')
plt.title("Top 10 Features Affecting Military Expenditure")
plt.xlabel("Feature Importance")
plt.show()

# Show predictions of the Random Forest model
predictions = rf_model.predict(X_test)

# Create a DataFrame to compare actual vs predicted values
predictions_df = pd.DataFrame({
    'Actual': y_test,
    'Predicted': predictions
})

predictions_df.head()

# Plot actual vs predicted values
plt.figure(figsize=(10, 6))
plt.scatter(y_test, predictions, alpha=0.6, color='blue')
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--')
plt.xlabel('Actual Expenditure')
plt.ylabel('Predicted Expenditure')
plt.title('Actual vs Predicted Expenditure - Random Forest Model')
plt.show()

# Correlation Analysis
df_corr = data[['Expenditure_Constant_2022', 'total_score', 'avg_rank', 'total_score_rank', 'Expenditure-Share_of_GDP']]
plt.figure(figsize=(10, 6))
sns.heatmap(df_corr.corr(), annot=True, cmap='coolwarm', linewidths=0.5)
plt.title("Correlation Between Conflict Metrics and Military Expenditure")
plt.show()

# Define the parameter grid for KMeans
param_grid = {
    'init': ['k-means++', 'random'],  # Initialization methods
    'n_init': [10, 20, 30],           # Number of centroid initializations
    'max_iter': [100, 500, 1000]       # Maximum number of iterations
}

# Custom grid search loop for KMeans
best_score = -1
best_params = {}
scores = []

for init in param_grid['init']:
    for n_init in param_grid['n_init']:
        for max_iter in param_grid['max_iter']:
            # Train KMeans with the current parameter combination
            kmeans = KMeans(n_clusters=optimal_clusters, init=init, n_init=n_init, max_iter=max_iter, random_state=42)
            labels = kmeans.fit_predict(reduced_data)
            
            # Compute the silhouette score
            score = silhouette_score(reduced_data, labels)
            scores.append((init, n_init, max_iter, score))
            
            # Update the best parameters if the current score is better
            # if score > best_score:
            best_score = score
            best_params = {'init': init, 'n_init': n_init, 'max_iter': max_iter}

# Train the final KMeans model with the best parameters
kmeans_best = KMeans(n_clusters=optimal_clusters, **best_params, random_state=42)
final_labels = kmeans_best.fit_predict(reduced_data)

# Visualize the final clusters
plt.figure(figsize=(10, 6))
plt.scatter(reduced_data[:, 0], reduced_data[:, 1], c=final_labels, cmap='viridis', s=50, alpha=0.7)
plt.title(f'KMeans Clustering with Optimized Parameters ({optimal_clusters} Clusters)')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.colorbar(label='Cluster')
plt.grid(True)
plt.show()

# Display the best parameters and silhouette score
results_df = pd.DataFrame(scores, columns=['init', 'n_init', 'max_iter', 'silhouette_score'])
best_results = pd.DataFrame([{'Best Parameters': best_params, 'Best Silhouette Score': best_score}])

# Print best parameters and score
print("Best KMeans Hyperparameters and Silhouette Score:")
print(best_results)

# Display all results
print("\nFull Hyperparameter Search Results:")
print(results_df)

Best KMeans Hyperparameters and Silhouette Score:
                                     Best Parameters  Best Silhouette Score
0  {'init': 'random', 'n_init': 30, 'max_iter': 1...               0.971119

Full Hyperparameter Search Results:
         init  n_init  max_iter  silhouette_score
0   k-means++      10       100          0.991880
1   k-means++      10       500          0.991880
2   k-means++      10      1000          0.991880
3   k-means++      20       100          0.991880
4   k-means++      20       500          0.991880
5   k-means++      20      1000          0.991880
6   k-means++      30       100          0.991880
7   k-means++      30       500          0.991880
8   k-means++      30      1000          0.991880
9      random      10       100          0.971119
10     random      10       500          0.971119
11     random      10      1000          0.971119
12     random      20       100          0.971119
13     random      20       500          0.971119
14     random      20      1000          0.971119
15     random      30       100          0.971119
16     random      30       500          0.971119
17     random      30      1000          0.971119

# 1. How does having a high conflict score affect national, regional, and subregional defense spending?
# Extracting the original categorical columns from the one-hot encoded dataframe
df = pd.get_dummies(data, columns=['Region', 'Subregion'], drop_first=True)
region_columns = [col for col in df.columns if 'Region_' in col]
subregion_columns = [col for col in df.columns if 'Subregion_' in col]

df['Original_Region'] = df[region_columns].idxmax(axis=1).str.replace('Region_', '')
df['Original_Subregion'] = df[subregion_columns].idxmax(axis=1).str.replace('Subregion_', '')

# Group by Original Region and Subregion and calculate average expenditure
region_exp = df.groupby('Original_Region')['Expenditure_Constant_2022'].mean()
subregion_exp = df.groupby('Original_Subregion')['Expenditure_Constant_2022'].mean()

print("Average Military Expenditure by Region (in Millions of Dollars):")
print(region_exp)

print("\nAverage Military Expenditure by Subregion (in Millions of Dollars):")
print(subregion_exp)

Average Military Expenditure by Region (in Millions of Dollars):
Original_Region
Americas          13180.149421
Asia & Oceania    13760.557746
Europe             8417.214983
Middle East       11057.362386
Name: Expenditure_Constant_2022, dtype: float64

Average Military Expenditure by Subregion (in Millions of Dollars):
Original_Subregion
Central Europe          1195.660545
East Asia              50240.653760
Eastern Europe         10805.445612
Middle East            11057.362386
North Africa            3384.418351
North America         416188.457695
Oceania                 6487.733853
South America           3918.272636
South Asia             11770.718225
South East Asia         2705.938999
Western Europe         13165.204485
sub-Saharan Africa       524.582137
Name: Expenditure_Constant_2022, dtype: float64

# 2. How does a country's GDP impact its defense spending?
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df, x='Expenditure-Share_of_GDP', y='Expenditure_Constant_2022')
plt.title("Impact of GDP Share on Military Expenditure")
plt.xlabel("Share of GDP in Military Expenditure")
plt.ylabel("Military Expenditure (Constant 2022)")
plt.show()

# 3. Does conflict score increase a country's spending as a greater proportion of GDP?
sns.lmplot(data=df, x='total_score', y='Expenditure-Share_of_GDP', hue='Original_Region', aspect=1.5)
plt.title("Impact of Conflict Score on Military Expenditure as % of GDP")
plt.xlabel("Total Conflict Score")
plt.ylabel("Military Expenditure as % of GDP")
plt.show()

# Visualizing Conflict Impact Across Subregions
plt.figure(figsize=(14, 8))
sns.boxplot(data=df, x='Original_Subregion', y='total_score', hue='Original_Subregion', palette='Set2', dodge=False)
plt.xticks(rotation=90)
plt.title("Conflict Score Across Subregions")
plt.xlabel("Subregion")
plt.ylabel("Total Conflict Score")
plt.legend().remove()
plt.show()

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

# Creating Heatmap for Conflict Scores Across Subregions
conflict_pivot = df.pivot_table(values='total_score', index='Original_Subregion', aggfunc='mean')
plt.figure(figsize=(12, 6))
sns.heatmap(conflict_pivot, annot=True, cmap='YlGnBu', linewidths=0.5)
plt.title("Heatmap of Average Conflict Scores Across Subregions")
plt.xlabel("Subregion")
plt.ylabel("Conflict Score")
plt.show()

# Creating Heatmap for GDP Share vs Conflict 
gdp_conflict_pivot = df.pivot_table(values='Expenditure-Share_of_GDP', index='Original_Subregion', columns='total_score', aggfunc='mean')
plt.figure(figsize=(16, 8))
sns.heatmap(gdp_conflict_pivot, annot=True, cmap='Reds', linewidths=0.5, fmt='.3f', cbar_kws={'label': 'GDP Share'})
plt.title("Heatmap of GDP Share vs Conflict Score Across Subregions")
plt.xlabel("Total Conflict Score")
plt.ylabel("Subregion")
plt.xticks(rotation=45)
plt.yticks(rotation=0)
plt.tight_layout()
plt.show()

# Analyzing Correlations Between GDP and Conflict
gdp_conflict_corr = df[['Expenditure-Share_of_GDP', 'total_score', 'avg_rank', 'total_score_rank']]
plt.figure(figsize=(10, 6))
sns.heatmap(gdp_conflict_corr.corr(), annot=True, cmap='viridis', linewidths=0.5)
plt.title("Correlation Between GDP Share and Conflict Metrics")
plt.show()

# Visualizing Defense Spending Trends Across Regions
plt.figure(figsize=(14, 8))
sns.lineplot(data=df, x='Year', y='Expenditure_Constant_2022', hue='Original_Region', marker='o', palette='Set1')
plt.title("Defense Spending Trends Across Regions")
plt.xlabel("Year")
plt.ylabel("Military Expenditure (Constant 2022)")
plt.xticks(rotation=45)
plt.legend(title='Region')
plt.tight_layout()
plt.show()

	Year	Expenditure-Share_of_Govt_spending	Expenditure-Share_of_GDP	Expenditure_Per_Capita	Expenditure_Constant_2022	index_level	total_rank_from_avg_rank	avg_rank	Deadliness_rank	Diffusion_rank	...	total_score_rank	Deadliness_raw	Diffusion_raw	Danger_raw	Fragmentation_raw	Deadliness_scaled	Diffusion_scaled	Danger_scaled	Fragmentation_scaled	total_score
count	3716.000000	3533.000000	3602.000000	3604.000000	3673.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	...	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000	3716.000000
mean	2011.556512	0.071187	0.019767	243.374084	11760.889052	3.491119	86.736006	71.252758	79.672497	43.389397	...	80.705867	918.312971	0.005629	252.864101	17.972282	0.022999	0.008146	0.039034	0.011832	0.081966
std	6.889732	0.060611	0.018137	426.197550	68296.745904	0.894834	49.421439	33.403444	40.775882	13.511428	...	42.837423	3629.622878	0.023324	710.539730	107.454813	0.090902	0.033754	0.109685	0.070740	0.231303
min	2000.000000	0.000672	0.000163	0.071803	0.000000	1.000000	2.000000	4.250000	2.000000	2.000000	...	2.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	2006.000000	0.031740	0.010229	17.077744	144.314628	3.000000	45.000000	42.750000	43.000000	37.000000	...	45.000000	0.000000	0.000000	2.000000	0.000000	0.000000	0.000000	0.000309	0.000000	0.001000
50%	2012.000000	0.051833	0.014996	69.809759	739.329049	4.000000	85.000000	76.500000	84.000000	51.000000	...	81.000000	7.000000	0.000000	13.000000	2.000000	0.000175	0.000000	0.002007	0.001317	0.005000
75%	2018.000000	0.091994	0.023754	270.439168	4523.330790	4.000000	131.000000	101.000000	125.000000	51.000000	...	117.000000	246.000000	0.001000	152.000000	9.000000	0.006161	0.001447	0.023464	0.005925	0.037000
max	2023.000000	0.581707	0.366531	5718.771025	990485.412100	4.000000	160.000000	114.750000	125.000000	51.000000	...	136.000000	37014.000000	0.176000	5426.000000	1519.000000	0.926995	0.254703	0.837604	1.000000	1.940000

	Country	Year	Expenditure_Constant_2022
619	Afghanistan	2004	183.648324
620	Afghanistan	2005	165.378054
621	Afghanistan	2006	167.189480
622	Afghanistan	2007	257.366345
623	Afghanistan	2008	224.292155

	Country	Linear R^2 Score	Random Forest R^2 Score
0	United States of America	-25.457937	-27.102278
1	China	0.694556	-5.880555
2	Russia	0.290156	-0.794920
3	Somalia	1.000000	0.000000
4	India	-0.237510	-10.619484
5	Saudi Arabia	-78.910826	-8.238869
6	United Kingdom	-0.080564	-3.584137
7	Ukraine	-0.674017	-0.720477
8	Germany	-11.194959	-7.560721
9	France	-1.866650	-2.222828
10	Japan	-2.388573	-1.605633
11	Korea, South	-3.369126	-19.531874
12	Italy	-10.028473	-2.944354
13	Australia	0.406479	-4.445582
14	Israel	-0.070533	-1.030049

	Country	Mean_Intercept	Feature_Importance	R_squared
0	United States of America	797284.510632	[1.0]	-94132.992031
1	China	134350.906012	[1.0]	-2626.882945
2	Russia	56886.987055	[1.0]	-277.411792
3	Somalia	126473.354300	[0.0]	-1786.684212
4	India	49049.911387	[1.0]	-163.153789
5	Saudi Arabia	57943.995012	[1.0]	-289.098608
6	United Kingdom	61806.560849	[1.0]	-296.795409
7	Ukraine	3475.918031	[1.0]	-29.979516
8	Germany	42074.764294	[1.0]	-89.737523
9	France	48054.356989	[1.0]	-138.682528
10	Japan	39764.717082	[1.0]	-73.304696
11	Korea, South	29782.014937	[1.0]	-22.886129
12	Italy	30678.855325	[1.0]	-29.425584
13	Australia	21684.055534	[1.0]	-1.814031
14	Israel	17716.323664	[1.0]	1.000000

	GovtSpending	GDP	PerCapita	ExpenditureConstant	TotalScore	Deadliness
0	-0.310629	-0.129417	-0.412730	0.059057	3.325245	1.587069
1	-0.389394	-0.005680	-0.420261	0.095522	3.325245	1.587069
2	-0.491608	-0.037074	-0.438759	0.099084	3.325245	1.587069
3	-0.591814	-0.256563	-0.456810	0.044395	3.325245	1.587069
4	-0.585301	-0.280179	-0.440198	0.052818	3.325245	1.587069

	Country	Year	Expenditure-Share_of_Govt_spending	Expenditure-Share_of_GDP	Expenditure_Per_Capita	Expenditure_Constant_2022	Region	Subregion	index_level	total_rank_from_avg_rank	...	Deadliness_scaled_rank	Diffusion_scaled_rank	Danger_scaled_rank	Fragmentation_scaled_rank	total_score_rank	Deadliness_raw	Diffusion_raw	Danger_raw	Fragmentation_raw	total_score
0	Brazil	2000	NaN	0.017307	64.500907	15643.43168	Americas	South America	1	6	...	10	31	3	9	6	6678	0.003	4117	67	0.851
1	Brazil	2001	0.047167	0.019519	61.332747	18119.43481	Americas	South America	1	6	...	10	31	3	9	6	6678	0.003	4117	67	0.851
2	Brazil	2002	0.041112	0.018958	53.550103	18361.33076	Americas	South America	1	6	...	10	31	3	9	6	6678	0.003	4117	67	0.851
3	Brazil	2003	0.035175	0.015035	45.956041	14647.82954	Americas	South America	1	6	...	10	31	3	9	6	6678	0.003	4117	67	0.851
4	Brazil	2004	0.035561	0.014613	52.945029	15219.74744	Americas	South America	1	6	...	10	31	3	9	6	6678	0.003	4117	67	0.851

	Country	Gradient Boosting R^2 Score
0	United States of America	-23.091954
1	China	-4.215515
2	Russia	-1.005995
3	Somalia	1.000000
4	India	-7.432400
5	Saudi Arabia	-10.015545
6	United Kingdom	-3.439904
7	Ukraine	-0.697479
8	Germany	-6.679121
9	France	-2.599453
10	Japan	-1.182671
11	Korea, South	-15.473089
12	Italy	-2.677620
13	Australia	-4.756625
14	Israel	-0.839813

	Actual	Predicted
3608	46897.50895	46986.231823
1140	20627.07245	23242.458146
3678	615755.93240	625818.797038
2935	69196.88482	75655.604037
3593	29352.65461	29520.882210

Forecasting Military Expenditure and Analyzing Spending Patterns Across Nations¶

1. Introduction:¶

2. Goal:¶

3. Dataset:¶

3.1 Package Requirements:¶

3.2 Imports:¶

3.3 Data Import:¶

4. Data Cleaning:¶

4.1 Imputing missing values¶

5. Exploratory Data Analysis:¶

Key Insights:¶

Key Insights:¶

6. Model Selection:¶

1. Supervised Learning¶

6.1.1 Linear Regression¶

6.1.2 Random Forest Regression¶

6.1.3 Gradient Boosting Regression¶

6.1.4 Prophet¶

2. Unsupervised Learning¶

Steps¶

1. Feature Selection:¶

2. Normalization:¶

3. Dimensionality Reduction:¶

4. Optimal Clusters Identification using Elbow Method Analysis:¶

6.2.1. K-Means Clustering:¶

6.2.2. Gaussian Mixture Model (GMM)¶

6.2.3. Agglomerative Hierarchical Clustering¶

6.2.4. DBSCAN (Density-Based Spatial Clustering of Applications with Noise)¶

Comparative Analysis:¶

Insights¶

Final Recommendation on Best Model Performance¶

7. Model Evaluation:¶

7.1 Hyper Parameter Tuning for Supervised Learning: Random Forest Regression¶

Summary¶

Summary:¶

Key Takeaways:¶

7.2 Hyper Parameter Tuning for Unsupervised Learning: K-Means¶

8. Model Analysis:¶

8.1. Supervised Learning using Random Forest Regression¶

Summary & Reasoning:¶

Reasoning:¶

Reasoning:¶

Reasoning:¶

Key Insights:¶

Reasoning:¶

8.2 Unsupervised Learning using K-Means Clustering¶

9. Conclusion:¶

10. Recommendations:¶

11. References:¶